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I. INTRODUCTION 



In general numerical simulations of the dynamics of soliton solutions in particle physics 
are of fundamental interest to test their stability. Specifically the dynamics of Einstein- Yang- 
Mills systems help to understand the interplay of gravity and non-abelian gauge fields. As 
analytical solutions for the fully nonlinear equations are usually not available in these cases, 
only numerical simulations can provide further insight. 

Bartnik and McKinnon |l| found that particle-like solutions for the Einstein- Yang-Mills 
system exist, while that is not the case for either pure Yang- Mills or the pure Einstein systems 
2J. Volkov et al. [3] investigated the Einstein- Yang-Mills equations with cosmological 
constant A to find an analogous class of solutions. 

In {4] Brodbeck et al. performed the respective stability analysis of these static Bartnick- 
McKinnon analogues. First steps in simulating the dynamical evolution of the the ground 
state n — 1 in the case of A = A reg (n) has been taken in jlj]. Different to the approach of [s| 
we use a Finite Element scheme for our numerical simulations and a reformulated differential 
equation. Special emphasis is on the boundary conditions in relation of the conservation of 
energy. In our approach we consider the static de Sitter manifold only. 

The paper is organized as follows. In section [Til we discuss the physical field equation that 
governs the dynamics of the monopole. In section III II we derive the initial- and boundary 
conditions in relation to the particle energy. The discretisation and related stability issues 
are discussed in section IIVI Finally numerical results are given in section |Vj 



II. THE PHYSICAL FIELD EQUATIONS 

In their analysis of Einstein- Yang-Mills fields with a cosmological constant Volkov et 
al. [3j] found that three different classes of solutions exist. These are static, spherically 
symmetric and are characterized by their value of A and the number of nodes n. The three 
different classes can be distinguished by their values of A. For A < A crit {n) the solutions 
approach asymptotically the de Sitter geometry. The solutions with A crit (n) < A < A reg (n) 
are the so called "bag of gold solutions". They have finite size, but still an horizon. All 
the solutions with non-vanishing A < A reg have non-vanishing magnetic charge. At A reg 
the solutions are compact solutions on a de Sitter manifold, where magnetic charges are no 
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longer defined. 

The smooth, spherically symmetric Yang-Mills field has a SU(2) gauge group. The com- 
ponents of the Gauge field in spherical coordinates for the Abelian gauge are 

A t =0, A#=--T 2 w(r,t), (la) 
A r =0, A v =- ^(-Tiw(r,t)sin^ + r 3 cos^), (lb) 

where g is the coupling constant of the Yang-Mills (YM) field and the Tj = ~<7j are the 
generators of the SU(2) group, <7j being the Pauli matrices. This gauge has been chosen 
as the Yang-Mills field in this case only depends on r and t, which reflects the spherical 
symmetry of the field. 



The de Sitter manifold is defined by its metric. In conformal coordinates 

by 

1 1 sin 2 r sin 2 r sin 2 1? 

sin t sin t sin t sin t 



it is given 



/ ± ± Olll / Dill / Olll U \ 

diag(^) = ( —2-, — t^t, — t^t, — ) • (2) 



The coordinates ranges are t G]0,7r[, r G [0,7r[, 1? G [0, 7r[ and <f G [0, 27r[. Its curvature 
is given by A/3, where A is the cosmological constant. The field equations can be derived 
through using the variational principle on the appropriate Lagrangian. The Yang-Mills 
Lagrangian is given by 

£> YM =-Tr(F tl „F^), with (3) 
=d ll A v -d v A ll + ig[A li ,A v \. (4) 

Fn U is the associated 2-form field of the gauge field and is also called the field strength tensor 



and [•, •] is the product of the su(2) ' 
Mills system have been derived in 
field on a static de Sitter manifold 



jie algebra. The equations for the coupled Einstein- Yang- 
5] , as well as the equations for the decoupled Yang-Mills 



,. w(w 2 - 1) n . . 

w - w " + — — 2 " = °- ( 5 ) 

sin r 

Here the dot denotes the derivative with respect to time and prime the derivative with 
respect to r. This is a semi-linear partial differential equation of hyperbolic type. 



3 



III. INITIAL AND BOUNDARY CONDITIONS 



In this paper we investigate the numerical solution of equation (J5]). It is defined on the 
compact domain 

Q = IxI, J=[0,7r] (R , 

and requires initial conditions for both w and w. For the initial conditions for w we will 
chose the stationary solution given in equation for the initial conditions for w we in- 
troduce a perturbation of the stationary system. Boundary conditions follow from energy 
considerations, see below. 

A. Initial Conditions 

As we have already mentioned before we are analysing the dynamical evolution of the 
static Bartnik-McKinnon analogue with n = 1 and A = A reg = 3/2. It is given by 

Wo = cos(r), r E I. (6) 

Thus it represents our initial condition for w. In the following we set the coupling constant 
g = 1. As the initial condition for w we choose 

w = sin(r)/oM) G L 2 (I), (7) 

with an arbitrary function fo, which was chosen in different numerical experiments differ- 
ently. 

B. Boundary Conditions 

The boundary conditions for the system will be derived from the finiteness of the energy. 
For systems on a Minkowski space-times, these types of boundary conditions are related 
to the asymptotic behavior of the field at infinity. These conditions are implemented by 
an appropriate choice of the functions-space, e.g. L 2 (uj), etc. As for the Minkowski space- 
times, for the compact de Sitter manifold the energy finiteness of w will be related to an 
appropriate choice of the function space. 
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1. Energy conditions for #, cp. 



As we are investigating a spherical symmetric Yang-Mills field the energy is independent 
of these variables. 



2. Energy Conditions for r. 



In [9j Chrusciel et al. derived a priori estimates of the solution space of Yang-Mills 
equations on globally hyperbolic four dimensional Lorentz manifolds. They show that the 
Yang-Mills gauge field belongs to the Sobolev space H 2 {R 3 ) with 

and an appropriate constant C(t) < oo. From the assumed symmetry (A^ is independent 
of if) and © it follows w(-,t) G H\I) n <C?' 1/2 (I) with 

\\w(;t)\\ mm <C(t). (8) 

The energy of the Yang-Mills field in the decoupled case is given by 

E(t)= [ T^CdV, 

J t=const. 

where = F^F u — ^g^F^F^ is the energy momentum tensor and n v = (sint, 0, 0, 0), 
= (1, 0, 0, 0) are the normal of the domain pointing outside and a Killing vector respec- 
tively. With the constant of proportionality C we rewrite the energy as 



E(t)=C j [w' 2 (r,t)+w 2 (r,t)^ 



w 2 {r,t)f 



2 sin 2 r 



dr 



l-w 2 (-,t) 



sin 



(9) 



L2. 



Due to the continuity of w(r, t) with respect to r, the finiteness of E(t) implies w(r, t) — ±1, 
Vr = 0, ii and t E I. We choose 



w(0,t) = 1, w(n,t) = -1, 



(10) 



Assuming the boundary conditions (JTUJ) and a field w(r,t) governed by equation dSJ), it is 
easily varified that the energy is conserved at all times. With ([6]) it follows 



E(t)=E(0) = C(3n/A+\\w \\l 2 ) 



(11) 
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As for the initial condition (J6]), we assume vanishing spatial derivatives of the solution at 
r = 0, 7r for all times. This follows from the assumption that the solution should describe 
a smooth solution on the three-dimensional sphere requiring tangential planes at the poles 
corresponding to r = 0, it. We therefore make the ansatz 

w(r, t) = cos(r) + sin(r)/(r, t), (r,t)eQ (12) 

with the unknown function f(r,t). From (jSJ) follows the equation for / 

sin(r)(/ - /") - 2 cos(r)/' + 2^^/ 

+3cos(r)/ 2 + sin(r)/ 3 = 0, (13) 

with the boundary and initial conditions 

f(0,t)=f(n,t)=0, Vtel, 

/(r,0) = 0,/(r,0) = / (r),V*eJ. (14) 

The following two terms in equation ( fl3l) need special attention: The scaling factor of 
d'Alembert operator (sin(r) vanishes at the boundary r = 0, 7r and the source term linear 
in / has infinite weights at the boundaries. The first difficulty is in view of a non-vanishing 
source at the boundary unproblematic. Assuming the same regularity for / as we have for 
w, i.e. f(-,t) E H\I) n C°' 1/2 (J), the linearly diver ging weight of the linear source term is 
also not problematic. 

IV. DISCRETISATION 

To integrate ( TTSl forwards in time, we transformed the equation to a set of first order 
equations, which we integrate by two different implicit schemes. The first is the classical 
implicit Euler methods, the second a modification of a two-stage Runge-Kutta method 



(Alexander scheme) [10|] , first introduced by in 11]. 

For the discretisation of the spatial dimension we chose the conforming Finite Element 
discretisation with continuous, piece-wise linear basis functions. We used a equidistant grid 
with mesh-size h = 7r/n r , n r > 0. The spatial integrals have been carried out with Maple 
with the exception of the fourth term in equation (fl3l) . This term cannot be integrated 
analytically and has been treated by appropriate integration rules (Simpson rule). The 



numerical approximation of these integrals has been investigated and chosen appropriately 
to avoid numerical artifacts. 

Additionally it was necessary to use Taylor expansions for the sine and cosine terms in 
the discretisation of the two non-linear terms in f. As an illustration for this issue we discuss 
the term cubic in f. Its discretisation is given by 

(sinr (fiipif,^), 

where (•, •) denotes the L 2 scalar product, ipi and if)j denote the piecewise linear basis func- 
tions. For j — it has the form Ci/q + C2/0 h + C3/0/1 2 + C4/1 • All the terms are similar, 
so we concentrate on 

C 4 = ^[-24cos(/i) + 6/isin(/i) + 24cos(2/i)-/i 3 sin(2/i) 
-Qh 2 cos(2/i) + 18/isin(2/i)]. 

Here special attention is required to control round-off errors. To accurately account for 
the difference 24 cos(2/i) — 24 cos(/i) an expansion of the trigonometric functions had to be 
carried out. If was found that an expansion up to the sixth order in h was sufficient, leading 
to an evaluation of C 4 according to 

^ h 2 67h 4 257/i 6 ^,, 8 , 

C A = 1 h 0(h s ). 

12 1680 43200 v ; 

Numerical experiments have shown that a further increase of the order were not necessary. 

V. NUMERICAL RESULTS AND CONCLUSIONS 

The simulations have been carried out with a software tool sEYMs-lD (simple Einstein- 



Yang-Mills solver for one dimension), a C++- code. A documentation is available with 12]. 

Several numerical experiments have been carried out. To validate to code, the conser- 
vation of energy has been considered. Its behavior over time is displayed in Fig. [T] on a 
sequence of successively refined grids. 

As a second numerical example we investigated the behavior of the field resulting from 
the initial conditions (fl4j) with 

/ = 10sin(r). (15) 
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Conserved quantity vs. time on different grids - space grid refinement 



I 
§- 

5 

c 

o 
U 



3.5345 



3.53445 



3.5344 - 



3.53435 - 



3.5343 



3.53425 



3.5342 



3.53415 



1 


1 1 


200x3200 






400x3200 






800x3200 






^\ 1600x3200 


- 




\ 3?OOy3?00 


1 


1 1 


1 1 1 



0.5 



1.5 



2.5 



3.5 



FIG. 1: Energy conservation of the YM field. The field f(r,t) has been computed using the 
Alexander scheme for time integration with initial conditions /o = sin(r). The label of the graphs 
indicate the spatial and temporal resolution n r x m, where n r and nt denote number of elements 
in the spatial and temporal direction, respectively. 



The amplitude has been chosen large enough in order to show the full non-linear dynamics 
of the system without destroying the original structure of the field. The simulations has 
been carried out on a 200 x 200 grid (using same notation as above) using the Implicit Euler 
method. It is displayed in Fig. [2j The result in terms of the original field w is shown in 
Fig. [3j The visualization has been carried out with the program xd3d 

As can be seen clearly, energy transfer takes place between the different modes resulting 
in a strongly perturbed steady-state solution. No plane wave solutions can be observed. 
This is a property intrinsic to Yang-Mills fields and has already been seen on flat manifolds. 
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FIG. 2: f(r,t) for initial conditions /o = 10sin(r). r, t and f(r,t) are varying along the x-, y- and 
z-axis, respectively 
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FIG. 3: w(r,t) for initial conditions /o = 10sin(r). r, t and w(r,t) are varying along the x-, y- and 
z-axis, respectively 
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